function [sag, cor, tra] = aft_SliceData3D(hfile, crd, opts)
% AFT::SliceData3D  - slice spatial data
%
% FORMAT:       [sag, cor, tra] = obj.SliceData3D(crd [, opts]);
%
% Input fields:
%
%       crd         1x3 coordinate to create oblique slices at
%       opts        optional settings
%        .conv      convention, either 'neurological' or {'radiological'}
%        .dir       direction, either of {'all'}, 'cor', 'sag', 'tra'
%        .fmr       required for AMR/MAP objects
%        .frame     min/max values (default: 256 box)
%        .mapvol    required for CMP/VMP objects, defaults to 1
%        .method    resampling method, one of {'cubic'}, 'linear', 'nearest'
%        .space     either of 'bvc', 'bvi', 'bvs', {'tal'}
%        .step      step size (default: 1)
%        .trans     4x4 transformation matrix, used for volume datatypes
%        .trf       TRF object(s) passed on to samplefmrspace
%        .v16       for VMR objects, if true access available V16 data
%
% Output fields:
%
%       sag         SAG slice data
%       cor         COR slice data
%       tra         TRA slice data
%
% TYPES: AMR, CMP, DMR, FMR, GLM, HDR, HEAD, MAP, MSK, SRF, VDW, VMP, VMR, VTC
%
% Note: for AMR/DMR/FMR/MAP objects, the trf option must be set with
%       at least the IA/FA (or FA alone) TRF object(s) to work

% Version:  v0.8a
% Build:    9102122
% Date:     Oct-21 2009, 10:08 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, true) || ...
   ~isa(crd, 'double') || ...
    numel(crd) ~= 3 || ...
    any(isinf(crd) | isnan(crd) | crd < -256 | crd > 256)
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'conv') || ...
   ~ischar(opts.conv) || ...
    isempty(opts.conv) || ...
   ~any(lower(opts.conv(1)) == 'nr')
    opts.conv = 'r';
else
    opts.conv = lower(opts.conv(1));
end
if ~isfield(opts, 'dir') || ...
   ~ischar(opts.dir) || ...
   ~any(strcmpi(opts.dir(:)', {'a', 'all', 'c', 'cor', 's', 'sag', 't', 'tra'}))
    opts.dir = 'a';
else
    opts.dir = lower(opts.dir(1));
end
if ~isfield(opts, 'fmr') || ...
    numel(opts.fmr) ~= 1 || ...
   (~isBVQXfile(opts.fmr, 'dmr') && ...
    ~isBVQXfile(opts.fmr, 'fmr'))
    opts.fmr = [];
end
if ~isfield(opts, 'frame') || ...
   ~isa(opts.frame, 'double') || ...
   ~isequal(size(opts.frame), [2, 3]) || ...
    any(isinf(opts.frame(:)) | isnan(opts.frame(:)) | opts.frame(:) < -256 | opts.frame(:) > 256)
    opts.frame = [];
end
if ~isfield(opts, 'mapvol') || ...
    numel(opts.mapvol) ~= 1 || ...
   ~isa(opts.mapvol, 'double') || ...
    isinf(opts.mapvol) || ...
    isnan(opts.mapvol) || ...
    opts.mapvol < 1
    opts.mapvol = 1;
end
if ~isfield(opts, 'method') || ...
   ~ischar(opts.method) || ...
   ~any(strcmpi(opts.method, {'cubic', 'linear', 'nearest'}))
    opts.method = 'cubic';
else
    opts.method = lower(opts.method(:)');
end
if ~isfield(opts, 'space') || ...
   ~ischar(opts.space) || ...
   ~any(strcmpi(opts.space(:)', {'bvc', 'bvi', 'bvs', 'tal'}))
    opts.space = 'tal';
else
    opts.space = lower(opts.space(:)');
end
if ~isfield(opts, 'step') || ...
   ~isa(opts.step, 'double') || ...
    isinf(opts.step) || ...
    isnan(opts.step) || ...
    opts.step <= 0 || ...
    opts.step >= 16
    opts.step = 1;
end
if ~isfield(opts, 'tal') || ...
    numel(opts.tal) ~= 1 || ...
   ~isBVQXfile(opts.tal, 'tal')
    opts.tal = {};
end
if ~isfield(opts, 'trans') || ...
   ~isa(opts.trans, 'double') || ...
   ~isequal(size(opts.trans), [4, 4]) || ...
    any(isinf(opts.trans(:)) | isnan(opts.trans(:))) || ...
    any(opts.trans(4, :) ~= [0, 0, 0, 1])
    opts.trans = [];
end
if ~isfield(opts, 'trf') || ...
   (~iscell(opts.trf) && ...
    (numel(opts.trf) ~= 1 || ...
     ~isBVQXfile(opts.trf, 'trf')))
    opts.trf = {};
elseif ~iscell(opts.trf)
    opts.trf = {opts.trf};
end
for tc = numel(opts.trf):-1:1
    if numel(opts.trf{tc}) ~= 1 || ...
       (~isBVQXfile(opts.trf{tc}, 'tal') && ...
        ~isBVQXfile(opts.trf{tc}, 'trf'))
        opts.trf(tc) = [];
    end
end
opts.trf = opts.trf(:)';
if ~isfield(opts, 'v16') || ...
   ~islogical(opts.v16) || ...
    numel(opts.v16) ~= 1
    opts.v16 = false;
end

% switch on filetype
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
ft = lower(sc.S.Extensions{1});

% for DMR/FMR/VTC objects, resolve transio, otherwise accept map
switch (ft)
    case {'amr'}
        if isempty(opts.fmr)
            error( ...
                'BVQXfile:MissingOption', ...
                'Option field .fmr must be set to access AMR in 3D space.' ...
            );
        end
        resetc = false;
        odata = zeros([size(bc.Slice(1).AMRData), numel(bc.Slice)]);
        for slc = 1:numel(bc.Slice)
            if istransio(bc.Slice(slc).AMRData)
                bc.Slice(slc).AMRData = resolve(bc.Slice(slc).AMRData);
                resetc = true;
            end
            odata(:, :, slc) = bc.Slice(slc).AMRData;
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
    case {'cmp'}
        nvol = numel(bc.Component);
        opts.mapvol = round(min(nvol, opts.mapvol));
        if istransio(bc.Component(opts.mapvol).CMPData)
            bc.Component(opts.mapvol).CMPData = ...
                resolve(bc.Component(opts.mapvol).CMPData);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Component(opts.mapvol).CMPData;
    case {'dmr'}
        if istransio(bc.DWIData)
            bc.DWIData = resolve(bc.DWIData);
            bvqxfile_setcont(hfile.L, bc);
        end
        switch bc.DataStorageFormat
            case {2}
                svol = 3;
            case {3}
                svol = 4;
            case {4}
                svol = 1;
            otherwise
                error( ...
                    'BVQXfile:IllegalSetting', ...
                    'Illegal DataStorageFormat setting in DMR object.' ...
                );
        end
        nvol = size(bc.DWIData, svol);
        opts.mapvol = round(min(nvol, opts.mapvol));
        switch (svol)
            case {1}
                odata = squeeze(bc.DWIData(opts.mapvol, :, :, :));
            case {3}
                odata = squeeze(bc.DWIData(:, :, opts.mapvol, :));
            case {4}
                odata = squeeze(bc.DWIData(:, :, :, opts.mapvol));
        end
    case {'fmr'}
        resetc = false;
        for slc = 1:numel(bc.Slice)
            if istransio(bc.Slice(slc).STCData)
                bc.Slice(slc).STCData = resolve(bc.Slice(slc).STCData);
                resetc = true;
            end
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
        szdt = size(bc.Slice(1).STCData);
        nvol = szdt(3);
        opts.mapvol = round(min(nvol, opts.mapvol));
        switch bc.DataStorageFormat
            case {1}
                odata = zeros([szdt(1:2), numel(bc.Slice)]);
                for slc = 1:numel(bc.Slice)
                    odata(:, :, slc) = bc.Slice(slc).STCData(:, :, opts.mapvol);
                end
            case {2}
                odata = squeeze(bc.Slice.STCData(:, :, opts.mapvol, :));
            otherwise
                error( ...
                    'BVQXfile:IllegalSetting', ...
                    'Illegal DataStorageFormat setting in DMR object.' ...
                );
        end
    case {'glm'}
        nvol = numel(bc.Predictor);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = aft_GetVolume(hfile, opts.mapvol);
    case {'hdr'}
        odata = bc.VoxelData;
        if istransio(odata)
            odata = resolve(odata);
            bc.VoxelData = odata;
            bvqxfile_setcont(hfile.L, bc);
        end
        
        % also get transformation matrix and set space to bvc!
        cfr = hdr_CoordinateFrame(hfile);
        if isempty(opts.trans)
            opts.trans = inv(cfr.Trf);
        else
            opts.trans = opts.trans * inv(cfr.Trf);
        end
        opts.frame = [128, 128, 128; -127.99, -127.99, -127.99];
        opts.space = 'bvc';
    case {'head'}
        nvol = numel(bc.Brick);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = bc.Brick(opts.mapvol).Data;
        if istransio(odata)
            odata = resolve(odata);
            bc.Brick(opts.mapvol).Data = odata;
            bvqxfile_setcont(hfile.L, bc);
        end
        
        % also get transformation matrix and set space to bvc!
        cfr = head_CoordinateFrame(hfile);
        if isempty(opts.trans)
            opts.trans = inv(cfr.Trf);
        else
            opts.trans = opts.trans * inv(cfr.Trf);
        end
        opts.frame = [128, 128, 128; -127.99, -127.99, -127.99];
        opts.space = 'bvc';
    case {'map'}
        if isempty(opts.fmr)
            error( ...
                'BVQXfile:MissingOption', ...
                'Option field .fmr must be set to access MAP in 3D space.' ...
            );
        end
        resetc = false;
        odata = zeros([size(bc.Map(1).Data), numel(bc.Map)]);
        for slc = 1:numel(bc.Map)
            if istransio(bc.Map(slc).Data)
                bc.Map(slc).Data = resolve(bc.Map(slc).Data);
                resetc = true;
            end
            odata(:, :, slc) = bc.Map(slc).Data;
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
    case {'msk'}
        if istransio(bc.Mask)
            bc.Mask = resolve(bc.Mask);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Mask;
    case {'srf'}
        if isempty(bc.VertexVMRData)
            srf_BackToVMR(hfile);
            bc = bvqxfile_getcont(hfile.L);
        end
        odata = bc.VertexVMRData;
    case {'vdw'}
        if istransio(bc.VDWData)
            bc.VDWData = resolve(bc.VDWData);
            bvqxfile_setcont(hfile.L, bc);
        end
        nvol = size(bc.VDWData, 1);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = squeeze(bc.VDWData(opts.mapvol, :, :, :));
    case {'vmp'}
        nvol = numel(bc.Map);
        opts.mapvol = round(min(nvol, opts.mapvol));
        if istransio(bc.Map(opts.mapvol).VMPData)
            bc.Map(opts.mapvol).VMPData = resolve(bc.Map(opts.mapvol).VMPData);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Map(opts.mapvol).VMPData;
        if bc.Map(opts.mapvol).EnableClusterCheck && ...
            isequal(size(odata), size(bc.Map(opts.mapvol).VMPDataCT))
            odata = odata .* bc.Map(opts.mapvol).VMPDataCT;
        end
    case {'vmr'}
        if ~isempty(bc.VMRData16) && ...
            opts.v16
            if istransio(bc.VMRData16)
                bc.VMRData16 = resolve(bc.VMRData16);
                bvqxfile_setcont(hfile.L, bc);
            end
            odata = bc.VMRData16;
        else
            if istransio(bc.VMRData)
                bc.VMRData = resolve(bc.VMRData);
                bvqxfile_setcont(hfile.L, bc);
            end
            odata = bc.VMRData;
        end
    case {'vtc'}
        if istransio(bc.VTCData)
            bc.VTCData = resolve(bc.VTCData);
            bvqxfile_setcont(hfile.L, bc);
        end
        nvol = size(bc.VTCData, 1);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = squeeze(bc.VTCData(opts.mapvol, :, :, :));
    otherwise
        error( ...
            'BVQXfile:ObjectTypeUnsupported', ...
            'Unsupported object type for ::SliceData3D method.' ...
        );
end

% check frame/stepsize
opts.step = opts.step([1, 1, 1]);
if isempty(opts.frame)
    switch (opts.space)
        case {'tal'}
            opts.frame = [128, 128, 128; -127.99, -127.99, -127.99];
        case {'bvi', 'bvs'}
            opts.frame = [0, 0, 0; 255.99, 255.99, 255.99];
        case {'bvc'}
            opts.frame = 1 + [0, 0, 0; (size(odata) - 0.01)];
    end
end
if opts.conv ~= 'r'
    opts.frame(:, 1) = opts.frame([2, 1], 1);
end
off = diff(opts.frame) < 0;
opts.step(off) = -opts.step(off);

% sample true BV spatial objects for nearest/linear in a fast way
if isempty(opts.trans) && ...
    any(strcmp(ft, {'cmp', 'glm', 'msk', 'vdw', 'vmp', 'vmr', 'vtc'}))
    
    % create three coord lists
    xc = opts.frame(1, 1):opts.step(1):opts.frame(2, 1);
    yc = opts.frame(1, 2):opts.step(2):opts.frame(2, 2);
    zc = opts.frame(1, 3):opts.step(3):opts.frame(2, 3);

    % get object layout
    if ~strcmp(ft, 'vmr')
        ol = [size(odata), 1, bc.XStart, bc.YStart, bc.ZStart, ...
            bc.XEnd, bc.YEnd, bc.ZEnd, bc.Resolution([1, 1, 1])];
    else
        ol = aft_Layout(hfile);
    end
    
    % depending on space
    switch (opts.space)
        case {'tal'}
            crd = (1 ./ ol([13, 11, 12])) .* ((128 - crd) - ol([7,5, 6]));
            xc = (1 / ol(13)) .* ((128 - xc) - ol(7));
            yc = (1 / ol(11)) .* ((128 - yc) - ol(5));
            zc = (1 / ol(12)) .* ((128 - zc) - ol(6));
            nstep = -opts.step(:)' ./ ol([13, 11, 12]);
        case {'bvs'}
            crd = (1 ./ ol([13, 11, 12])) .* (crd - ol([7, 5, 6]));
            xc = (1 / ol(13)) .* (xc - ol(7));
            yc = (1 / ol(11)) .* (yc - ol(5));
            zc = (1 / ol(12)) .* (zc - ol(6));
            nstep = opts.step(:)' ./ ol([13, 11, 12]);
        case {'bvi'}
            crd = (1 ./ ol([13, 11, 12])) .* crd([3, 1, 2]);
            xs = xc;
            xc = (1 / ol(13)) .* (zc - ol(7));
            zc = (1 / ol(12)) .* (yc - ol(6));
            yc = (1 / ol(11)) .* (xs - ol(5));
            nstep = opts.step([3, 1, 2]) ./ ol([13, 11, 12]);
        case {'bvc'}
            crd = crd([3, 1, 2]);
            xs = xc;
            xc = zc;
            zc = yc;
            yc = xs;
            nstep = opts.step([3, 1, 2]);
    end
    
    % override method if possible
    if all(abs(nstep) == 1) && ...
        xc(1) == round(xc(1)) && ...
        yc(1) == round(yc(1)) && ...
        zc(1) == round(zc(1))
        opts.method = 'nearest';
    end
    
    % prepare outputs
    if opts.dir == 'a'
        sag = zeros(numel(zc), numel(yc));
        cor = zeros(numel(zc), numel(xc));
        tra = zeros(numel(yc), numel(xc));
    elseif opts.dir == 's'
        sag = zeros(numel(zc), numel(yc));
    elseif opts.dir == 'c'
        sag = zeros(numel(zc), numel(xc));
    else
        sag = zeros(numel(yc), numel(xc));
    end
    
    % sample only within valid coordinates
    xv = (xc > -0.5) & (xc < (size(odata, 3)) - 0.5);
    yv = (yc > -0.5) & (yc < (size(odata, 1)) - 0.5);
    zv = (zc > -0.5) & (zc < (size(odata, 2)) - 0.5);
    xf = xc(findfirst(xv));
    yf = yc(findfirst(yv));
    zf = zc(findfirst(zv));
    if isempty(xf) || ...
        isempty(yf) || ...
        isempty(zf)
        return;
    end
    
    % now sample only required parts, nearest
    if strcmpi(opts.method, 'nearest')
        crd = 1 + round(crd);
        xc = 1 + round(xc(xv));
        yc = 1 + round(yc(yv));
        zc = 1 + round(zc(zv));
        if crd(1) > 0 && ...
            crd(1) <= size(odata, 3) && ...
            any('as' == opts.dir)
            sag(zv, yv) = odata(yc, zc, crd(1))';
        end
        if crd(2) > 0 && ...
            crd(2) <= size(odata, 1)
            if opts.dir == 'a'
                cor(zv, xv) = odata(crd(2), zc, xc);
            elseif opts.dir == 'c'
                sag(zv, xv) = odata(crd(2), zc, xc);
            end
        end
        if crd(3) > 0 && ...
            crd(3) <= size(odata, 2)
            if opts.dir == 'a'
                tra(yv, xv) = odata(yc, crd(3), xc);
            elseif opts.dir == 't'
                sag(yv, xv) = odata(yc, crd(3), xc);
            end
        end
        
    % or highed interp
    else
        crd = crd + 1;
        xt = xf + (sum(xv) - 0.25) * nstep(1);
        yt = yf + (sum(yv) - 0.25) * nstep(2);
        zt = zf + (sum(zv) - 0.25) * nstep(3);
        
        crdf = [Inf, Inf, Inf; yf, zf, xf; nstep([2, 3, 1]); yt, zt, xt];
        crdf([2, 4], :) = 1 + crdf([2, 4], :);
        if crd(1) > -0.5 && ...
            crd(1) < size(odata, 3) && ...
            any('as' == opts.dir)
            crds = crdf;
            crds([2, 4], 3) = crd(1);
            sag(zv, yv) = flexinterpn_method(odata, crds, 0, opts.method)';
        end
        if crd(2) > -0.5 && ...
            crd(2) < size(odata, 1)
            crds = crdf;
            crds([2, 4], 1) = crd(2);
            if opts.dir == 'a'
                cor(zv, xv) = flexinterpn_method(odata, crds, 0, opts.method);
            elseif opts.dir == 'c'
                sag(zv, xv) = flexinterpn_method(odata, crds, 0, opts.method);
            end
        end
        if crd(3) > -0.5 && ...
            crd(3) < size(odata, 2)
            crds = crdf;
            crds([2, 4], 2) = crd(3);
            if opts.dir == 'a'
                tra(yv, xv) = flexinterpn_method(odata, crds, 0, opts.method);
            elseif opts.dir == 't'
                sag(yv, xv) = flexinterpn_method(odata, crds, 0, opts.method);
            end
        end
    end

% sample FMR based data differently
elseif any(strcmp(ft, {'amr', 'dmr', 'fmr', 'map'}))
    
    % build 1-d vectors, then pairs of coords
    xv = opts.frame(1, 1):opts.step:opts.frame(2, 1);
    yv = opts.frame(1, 2):opts.step:opts.frame(2, 2);
    zv = opts.frame(1, 3):opts.step:opts.frame(2, 3);
    [ys, xs] = ndgrid(zv, yv);
    [yc, xc] = ndgrid(zv, xv);
    [yt, xt] = ndgrid(yv, xv);
    if opts.dir == 'a'
        crd = [[crd(1) .* ones(numel(xs), 1), xs(:), ys(:)]; ...
               [xc(:), crd(2) .* ones(numel(xc), 1), yc(:)]; ...
               [xt(:), yt(:), crd(3) .* ones(numel(xt), 1)]];
    elseif opts.dir == 's'
        crd = [crd(1) .* ones(numel(xs), 1), xs(:), ys(:)];
    elseif opts.dir == 'c'
        crd = [xc(:), crd(2) .* ones(numel(xc), 1), yc(:)];
    else
        crd = [xt(:), yt(:), crd(3) .* ones(numel(xt), 1)];
    end
    
    % convert coordinates to BV system
    if strcmp(opts.space, 'tal')
        crd = 128 - crd;
    elseif strcmp(opts.space, 'bvi')
        crd = crd(:, [3, 1, 2]);
    elseif ~strcmp(opts.space, 'bvs')
        error( ...
            'BVQXfile:InvalidOption', ...
            'Sampling of DMR/FMR based data must not have bvc space.' ...
        );
    end
    
    % use AMR/MAP option DMR/FMR
    if ~isempty(opts.fmr)
        hfile = opts.fmr;
    end
    
    % sample data with samplefmrspace function
    data = samplefmrspace(odata, crd, hfile, ...
        [opts.trf, opts.tal], opts.method);
    
    % split into sag, cor, tra
    if any('as' == opts.dir)
        sag = reshape(data(1:numel(xs)), size(xs));
    end
    if opts.dir == 'a'
        cor = reshape(data(numel(xs)+1:numel(xs)+numel(xc)), size(xc));
        tra = reshape(data(numel(xs)+numel(xc)+1:end), size(xt));
    elseif opts.dir == 'c'
        sag = reshape(data(1:numel(xc)), size(xc));
    elseif opts.dir == 't'
        sag = reshape(data(1:numel(xt)), size(xt));
    end
    
% sample "already 3D" formats
else
    
    % get bounding box for valid objects and get transformation
    if ~strcmp(opts.space, 'bvc')
        ftrf = {bvcoordconv(zeros(0, 3), [opts.space '2bvc'], aft_BoundingBox(hfile))};
    else
        ftrf = {};
    end
    
    % any additional transformation?
    if ~isempty(opts.trans)
        if isempty(ftrf)
            ftrf = {opts.trans};
        else
            ftrf = {ftrf{1} * opts.trans};
        end
    end

    % sample data
    if any('as' == opts.dir)
        sf = opts.frame;
        sf(:, 1) = crd(1);
        sag = squeeze(flexinterpn_method(odata, ...
            [Inf, Inf, Inf; sf(1, :); opts.step; sf(2, :)], ...
            0, ftrf{:}, opts.method))';
    end
    sf = opts.frame;
    sf(:, 2) = crd(2);
    if opts.dir == 'a'
        cor = squeeze(flexinterpn_method(odata, ...
            [Inf, Inf, Inf; sf(1, :); opts.step; sf(2, :)], ...
            0, ftrf{:}, opts.method))';
    elseif opts.dir == 'c'
        sag = squeeze(flexinterpn_method(odata, ...
            [Inf, Inf, Inf; sf(1, :); opts.step; sf(2, :)], ...
            0, ftrf{:}, opts.method))';
    end
    sf = opts.frame;
    sf(:, 3) = crd(3);
    if opts.dir == 'a'
        tra = squeeze(flexinterpn_method(odata, ...
            [Inf, Inf, Inf; sf(1, :); opts.step; sf(2, :)], ...
            0, ftrf{:}, opts.method))';
    elseif opts.dir == 't'
        sag = squeeze(flexinterpn_method(odata, ...
            [Inf, Inf, Inf; sf(1, :); opts.step; sf(2, :)], ...
            0, ftrf{:}, opts.method))';
    end
end
